System and method for a hierarchical Bayesian-map approach for solving inverse problems

ABSTRACT

A method of reconstructing an image of an object, the method including: determining, by a plurality of sensors, a waveform based on the object, wherein the plurality of sensors view the object from a plurality of directions; determining, by a pre-processing module, a plurality of measurements of the object using the waveform, wherein the plurality of measurements are arranged in a vector form; determining, by an option module, a sampling matrix, a dictionary, and a noise factor, wherein the sampling matrix represents a geometric arrangement of the plurality of sensors, and the dictionary is pre-selected by the option module; estimating, by an estimation module, a coefficient vector using the measurements, the sampling matrix, and the noise factor; and reconstructing, by a reconstruction module, the image, using the coefficient vector and the dictionary.

GOVERNMENT INTEREST

The embodiments described herein may be manufactured, used, and/orlicensed by or for the United States Government without the payment ofroyalties thereon.

BACKGROUND

Technical Field

The embodiments herein relate to imaging, and more particularly toreconstructing images using one or more sensors.

Description of the Related Art

Imaging refers to a class of inverse problems wherein the centralobjective is to form the image of an object or scene of interest that istypically being sensed by one more sensors each furnishing complementarybut correlated sources of information, and each of which is potentiallycontaminated by noise and distortion. What distinguishes this from themore general class of inverse and regression problems is that imagesoriginating from typical empirical sources (such as natural images) havea special structure that makes it possible to inject informed priormodels into the inference process that can enhance the quality of thereconstructed image. Accordingly, a wide variety of analytical andstatistical approaches have proliferated over the past several decadesin diverse imaging applications such as radar, medical, and astronomicalimaging.

SUMMARY

In view of the foregoing, an embodiment herein provides a method ofreconstructing an image of an object, the method comprising:determining, by a plurality of sensors, a waveform based on the object,wherein the plurality of sensors view the object from a plurality ofdirections; determining, by a pre-processing module, a plurality ofmeasurements of the object using the waveform, wherein the plurality ofmeasurements are arranged in a vector form; determining, by an optionmodule, a sampling matrix, a dictionary, and a noise factor, wherein thesampling matrix represents a geometric arrangement of the plurality ofsensors, and the dictionary is pre-selected by the option module;estimating, by an estimation module, a coefficient vector using themeasurements, the sampling matrix, and the noise factor; andreconstructing, by a reconstruction module, the image, using thecoefficient vector and the dictionary.

The estimating the coefficient vector may comprise computing, by a firstvariable module: a first variable, using a pre-selected non-linearfactor; and a multi-scale Gaussian tree structure, using a quad-treedecomposition of the image, the sampling matrix, the dictionary, and themeasurements. Estimating the coefficient vectors may further comprise:estimating, by a parameter module, a structural parameter based on aparent-child relationship for each node in the tree structure;repeating, by a loop module: the computing of the first variable and themulti-scale Gaussian tree structure, and the estimating of thestructural parameter, across the tree structure, until the structuralparameter is lower than a first pre-selected threshold. The estimatingthe coefficient vectors may thrther comprise: computing, by a secondvariable module, a second variable based on the first variable, thesampling matrix, a variable selection operator, and a secondpre-selected threshold; and computing, by a coefficient module, thecoefficient vector based on a Hadamard product of the first variable andthe second variable.

Estimating the coefficient vector may comprise: initializing, by theestimation module, x⁰=|h⁻¹(Ψ^(T)y)| and n=0, wherein x⁰ is an initialapproximation of a temporary parameter, h is a nonlinear factor, Ψ isthe sampling matrix, y is a vector comprising the measurements, and n isa loop counter; calculating, by the estimation module, a descentdirection d^(n); determining, by the estimation module, a step size λ;and computing, by the estimation module, x^(n+1)=x^(n)+λd^(n), whereinx^(n) is an n^(th) approximation for the temporary variable and x^(n+1)is an (n+1)^(th) approximation for the temporary variable.

The method may further comprise: incrementing, by the estimation module,the loop counter n; computing, by the estimation module, x*, byrepeating the calculating the descent direction d^(n), the determiningthe step size λ, the computing x^(n+1), and the incrementing the loopcounter n until a norm of a steepest descent vector is smaller than apre-selected threshold, wherein x* is a temporary variable. The methodmay further comprise: computing, by the estimation module, z*=h(x*);calculating, by the estimation module,

=diag(S_(λ)[z*]) and Λ_(R)=diag (z*), where

$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix};} \right.}$calculating, by the estimation module, u*=L⁻¹R where L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is the noise factor; and calculating, by the estimation module, c*=

u* where c*ε

^(n) is the coefficient vector. Reconstructing the image may furthercomprise: calculating, by the reconstruction module, I=Φc* where Iε

^(n) is the reconstructed image and Φε

^(d×n) is the pre-selected dictionary. The plurality of sensors may beindependent from each other, wherein the dictionary may comprise atleast one class of wavelet dictionaries, and wherein the sampling matrixmay comprise a sampling operator determined by a transmitted waveformassociated with the measurements and the geometric arrangement of theplurality of sensors.

Another embodiment herein provides an imaging device comprising: aplurality of sensors configured to generate a waveform based on anobject; a pre-processor configured to determine a plurality ofmeasurements using the waveform, wherein the plurality of measurementsare arranged in a vector form; and a central imaging device comprising:an option module configured to determine a sampling matrix, apre-selected dictionary, a pre-selected non-linear factor, a first preselected threshold, a second pre-selected threshold, and a noise factor,wherein option module determines the sampling matrix using a geometricarrangement of the plurality of sensors; an estimation module configuredto estimate a coefficient vector using the plurality of measurements,the sampling matrix, and the pre-selected dictionary; and areconstruction module configured to reconstruct an image of the object,using the coefficient vector and the pre-selected dictionary.

The estimation module may further comprise: a first variable moduleconfigured to: compute a first variable using a pre-selected non-linearfactor; and determine a multi-scale Gaussian tree structure using aquad-tree decomposition of the electronic source image model, thesampling matrix, the pre-selected dictionary, and the measurements. Theestimation module may further comprise: a parameter module configured todetermine a structural parameter using a parent-child relationship foreach node in the multi-scale Gaussian tree structure; and a loop moduleconfigured to control operation of the first variable module and theparameter module across the multi-scale Gaussian tree structure untilthe structural parameter is lower than the first pre-selected threshold.The device may further comprise: a second variable module configured todetermine a second variable using the first variable, the samplingmatrix, the variable selection operator and the second pre-selectedthreshold; and a coefficient module configured to determine thecoefficient vector using a Hadamard product of the first variable andthe second variable.

The estimation module may be further configured to: initialize x⁰=|h³¹¹(Ψ^(T)y)| and n=0, wherein x⁰ is an initial approximation of atemporary parameter, h is a nonlinear factor, Ψ is the sampling matrix,y is a vector comprising the measurements, and n is a loop counter;calculate a descent direction d^(n); determine a step size λ; andcompute x^(n+1)=x^(n)+λd^(n), wherein x^(n) is an n^(th) approximationfor the temporary variable and x^(n+1) is an or (n+1)^(th) approximationfor the temporary variable. The estimation module may be furtherconfigured to: increment the loop counter n; compute x* by repeating thecalculating the descent direction d^(n), the determining the step sizeλ, the computing x^(n+), and the incrementing the loop counter n until anorm of a steepest descent vector is smaller than a pre-selectedthreshold. The estimation module may be further configured to: computez*=h(x*); calculate

=diag(S_(λ)[z*]) and Λ_(R)=diag(z*) where

$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix};} \right.}$calculate u*=L⁻¹R where L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is the noise factor; and calculate c*=z*⊙u* where c*ε

^(n) is the coefficient vector. The reconstruction module may be furtherconfigured to reconstruct the image by calculating I=Φc* where Iε

^(n) is the reconstructed image and Φε

^(d×n) is the pre-selected dictionary.

Another embodiment provides a non-transitory program storage devicereadable by computer, tangibly embodying a program of instructionsexecutable by the computer to perform a method for reconstructing animage of an object, the method comprising: detemining, by a plurality ofsensors, a waveform based on the object, wherein the plurality ofsensors view the object from a plurality of directions; determining, bya pre-processing module, a plurality of measurements of the object usingthe waveform, wherein the plurality of measurements are arranged in avector form; determining, by an option module, a sampling matrix, adictionary, and a noise factor, wherein the sampling matrix represents ageometric arrangement of the plurality of sensors, and the dictionary ispre-selected by the option module; estimating, by an estimation module,a coefficient vector using the measurements, the sampling matrix, andthe noise factor; and reconstructing, by a reconstruction module, theimage, using the coefficient vector and the dictionary.

Estimating the coefficient vector may comprise computing, by a firstvariable module a first variable, using a pre-selected non-linearfactor; and a multi-scale Gaussian tree structure, using a quad-treedecomposition of the image, the sampling matrix, the dictionary, and themeasurements. Estimating the coefficient vectors may further comprise:estimating, by a parameter module, a structural parameter based on aparent-child relationship for each node in the tree structure;repeating, by a loop module: the computing of the first variable and themulti-scale Gaussian tree structure, and the estimating of thestructural parameter, across the tree structure, until the structuralparameter is lower than a first pre-selected threshold. Reconstructingthe image may further comprise using global compound Gaussian model as astatistical prior. Reconstructing the image may further comprise usinghierarchical Bayesian maximum a posteriori and using global compoundGaussian model as a statistical prior.

These and other aspects of the embodiments herein will be betterappreciated and understood when considered in conjunction with thefollowing description and the accompanying drawings. It should beunderstood, however, that the following descriptions, while indicatingpreferred embodiments and numerous specific details thereof, are givenby way of illustration and not of limitation. Many changes andmodifications may be made within the scope of the embodiments hereinwithout departing from the spirit thereof, and the embodiments hereininclude all such modifications.

BRIEF DESCRIPTION OF THE DRAWINGS

The embodiments herein will be better understood from the followingdetailed description with reference to the drawings, in which:

FIG. 1 is a schematic diagram illustrating a an imaging system,according to an embodiment herein;

FIG. 2 is a schematic block diagram illustrating an imaging device,according to an embodiment herein;

FIG. 3 is a schematic block diagram illustrating an imaging device,according to an embodiment herein;

FIG. 4 is a flowchart illustrating a method, according to an embodimentherein;

FIG. 5 is a flowchart illustrating a method, according to an embodimentherein;

FIG. 6 is a flowchart illustrating a method, according to an embodimentherein;

FIG. 7 is a schematic diagram illustrating a quad-tree structure forcapturing a non-Gaussian interaction between wavelet coefficients,according to an embodiment herein;

FIG. 8 is a graph illustrating simulation results, according to anembodiment herein;

FIG. 9 is a graph illustrating simulation results, according to anembodiment herein;

FIG. 10 is set of photographic images, including an original image andseveral images reconstructed according to an embodiment herein;

FIG. 11 is a set of photographic images, including an original image andseveral images reconstructed according to an embodiment herein; and

FIG. 12 is a schematic diagram illustrating an exemplary computerarchitecture used in accordance with the embodiments herein.

DETAILED DESCRIPTION

The embodiments herein and the various features and advantageous detailsthereof are explained more fully with reference to the non-limitingembodiments that are illustrated in the accompanying drawings anddetailed in the following description. Descriptions of well-knowncomponents and processing techniques are omitted so as to notunnecessarily obscure the embodiments herein. The examples used hereinare intended merely to facilitate an understanding of ways in which theembodiments herein may be practiced and to further enable those of skillin the art to practice the embodiments herein. Accordingly, the examplesshould not be construed as limiting the scope of the embodiments herein.

The embodiments herein provide a Hierarchical Bayesian-MAP (Maximum aPosteriori) method for solving inverse problems. An embodiment hereinreconstructs an image that is viewed from multiple directions byindependent sensors, subject to a global Compound Gaussian (CG) priorplaced on the signal to be estimated. An embodiment herein approachesthis inverse problem in imaging based on a Hierarchical Bayesian-MAP(HB-MAP) formulation. An embodiment herein describes the basic inverseproblem of multi-sensor (tomographic) imaging. Given the measurementsrecorded by the sensors, a problem may be reconstructing the image witha high degree of fidelity. An embodiment uses a Probabilistic GraphicalModeling extension of the CG distribution as a global image prior into aHierarchical Bayesian inference procedure. In an embodiment, the globalCG model incorporated within the HB-MAP inference procedure is used forsolving inverse problems and imaging. The HB-MAP algorithm solves adifficult non-convex optimization problem. The Monte-Carlo simulations,described herein, show that that this approach works for over a broadrange of conditions—high and low sparsity cases—whereas ComprehensiveSensing (CS) and sparsity based approaches work only for high-sparsitycases. An embodiment herein provides a way to utilize the global CGmodel which subsumes the Laplacian—the basis of compressive sensing andmany other statistical models—for solving, difficult inverse problemssuch as imaging. An embodiment herein represents anumerical-technological advancement to realize the application of globalCG models for imaging in practice. Embodiments herein can be used forradar imaging (SAR, ISAR etc.), medical imaging, and a broad class ofinverse problems involving empirically derived data sources.

Referring now to the drawings, and more particularly to FIGS. 1 through12, where similar reference characters denote corresponding featuresconsistently throughout the figures, there are shown preferredembodiments.

FIG. 1 is a schematic diagram illustrating an imaging system 100according to an embodiment herein. Imaging system 100 may includeimaging device 150 that receives the signals resulting from thereflection of waveforms 152 (transmitted and received at diverse sensingangles) from object/scene of interest 154. Imaging system 100 can alsoinclude pre-processor 156 that conditions the signals by factoring outsignal variations that are due to various systemic factors.Pre-processor 156 produces data that is processed by central imagingdevice 158 to generate estimated image 160. Pre-processor 156 canencapsulate much of the system- and domain-level expertise associatedwith the application domain of interest. For example, in radar imaging,typical pre-processing steps can include pulse compression (generalizedmatched filtering) and range alignment (which corrects for the differenttime delays associated with the transmission of pulse sequences throughthe propagation media) which collectively ensure that the effectivesystem can be approximately modeled as a linear shift-invarianttransformation-thus rendering the resulting signals amenable toefficient Fourier analysis tools.

An embodiment herein provides optimum image reconstruction conditionedon the aforementioned pre-processing operations. An effective forwardmodel can be captured by the following linear system model:y=ΨΦc+v  (1)where:

-   yε    ^(m) is the measurements in vector form,-   Ψε    ^(m×d) is the effective measurement matrix,-   I=Φcε    ^(n) is the vectorized image of interest,-   Φε    ^(d×n) is the dictionary in which one may choose to represent the    unknown image I,-   cε    ^(n) is the unknown coefficient vector which is the object of the    inference methodology,    and-   vε    ^(m) is the effective measurement noise.

An embodiment herein solves equation (1) by estimating the coefficientvector c given the observed measurements y. The dictionary Φ, withrespect to which the image is represented, is chosen a priori and canbe, for example, but not limited to, any class of wavelet dictionaries.The dictionary can also be a member of a more general class ofover-complete dictionaries. Embodiments herein can be applied to anytype of sampling operator Ψ, for example, but not limited to, the Radontransform which underlies some imaging applications. In practice, thematrix Ψ may be determined by both the transmitted waveform and thegeometric arrangement of the sensors in space and time.

Different approaches may be used to solve for c in Eq. (1). For imagingapplications, the traditional approach is the Filtered Back-projection(FBP) methodology (Deans) wherein the underlying dictionary Φ includesFourier atoms. A variety of Compressive-Sensing(CS)/Sparsity-based-reconstruction approaches including Large-ScaleIl-Regularized Least Squares (Il-Is), Orthogonal Matching Pursuit (OMP),Regularized Orthogonal Matching Pursuit (ROMP), Iterative signalrecovery from incomplete and inaccurate samples (CoSaMP), and BayesianCompressive Sensing (BCS) have been applied for this purpose. Anembodiment herein encompasses CS and sparse reconstruction approachesand treat them as special cases.

From Eq. (1) an embodiment herein provides a Bayesian-Maximum aposteriori (MAP) estimate of c. The prefix “Bayesian” emphasizes thatprobabilities in embodiments herein are assigned without a correspondingnotion of sampling or frequency. Furthermore an embodiment hereinproduces covariance estimates which could be used to determineconfidence intervals.

The Bayesian-MAP estimate of c is given as follows:

$\begin{matrix}{c^{*} = {{{argmax}_{c}\log\;{P\left( {c\text{|}y} \right)}} \propto}} & {{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}(2)} & {{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}(2)} \\{{argmax}_{c}\left\{ {{\log\;{P\left( {y\text{|}c} \right)}} + {\log\;{P(c)}}} \right\}} & {(3)} & {(3)} \\{= {{argmin}_{c}\left\{ {{{y - {{\Psi\Phi}c}}}_{2}^{2} - {\log\;{P(c)}}} \right\}}} & {(4)} & {(4)}\end{matrix}$where, in Eq. (4) v˜

(0,Σ_(v)=1).

The distribution P(c) encapsulates the statistical prior knowledge aboutthe scene structure. For the specific choice of a Laplacian (i.e.P(c)=exp(−λ∥c∥₁)), Eq. (2) reduces to the LASSO methodology or the BasisPursuit methodology:c*=argmin _(c) {∥y−ΨΦc∥2/2+λ∥c∥ ₁}  (5)

The Laplacian distribution is not sufficiently rich for modeling thestatistics of wavelet coefficients of natural scenes when sensed byoptical or radar sensors. In an embodiment herein, a probabilisticgraphical modeling extension of the Compound Gaussian (CG) distributionis used as a candidate for P(c). This extension subsumes distributionssuch as, for example, but not limited to, the Laplacian, the generalizedGaussian, and the alpha-stable distribution.

The Bayesian-MAP estimation problem, under the CG modeling, lends itselfto a Hierarchical Bayes formulation that can yield superior performanceto traditional CS methodologies. In Bayesian statistics, a maximum aposteriori probability (MAP) estimate is a mode of the posteriordistribution. The MAP can be used to obtain a point estimate of anunobserved quantity on the basis of empirical data. It is closelyrelated to Fisher's method of maximum likelihood (ML), but employs anaugmented optimization objective which incorporates a prior distributionover the quantity one wants to estimate.

In Eq. (1), c can be modeled as a random vector that can be decomposedinto the following Hadamard product form:c=z⊙u  (6)such that:i) u˜

(0,P_(u))

-   -   z=h(x)        x follows a Multi-scale Gaussian Tree structure        ii) u and z are independent random variables        iii)        [z²]=1        iv) h is a non-linearity (which ultimately controls the sparse        structure of c).

FIG. 2, with reference to FIG. 1, is a schematic diagram illustratingthe imaging device 150 for reconstructing an electronic source image ofthe target that is viewed by independent sensors, according to anembodiment herein. In an embodiment, imaging device 150 includes sensors184 configured to generate the waveforms 152 based on an object 180. Inan embodiment, imaging device 150 includes the pre-processor 156 and thecentral imaging device 158. Pre-processor 156 processes thebackscattered waveforms 190 and furnishes the resulting measurements tothe central imaging device 150. The pre-processor 156 is furtherconfigured to determine measurements 174 of the object/scene beingsensed. Measurements 174 may be in a vector form.

The measurements 174 may be in a vector form. Imaging device 150 mayalso include options module 165 configured to determine a samplingmatrix 173, a pre-selected dictionary 161, and a noise factor 166.Sampling matrix 173 can be based on a geometric arrangement of sensors184 collecting data of an object 180. Imaging device 150 may furtherinclude estimation module 167 configured to estimate coefficient vector162 based on measurements 174, sampling matrix 173, pre-selecteddictionary 161, and noise factor 166.

Imaging device 150 may further include reconstruction module 169configured to prepare, a reconstructed image 160 from the measurements174. In an embodiment, reconstruction module 169 uses any of dictionary161 and estimated coefficient vector 162 to generate the reconstructedimage 160. In an embodiment, dictionary 161 may include at least oneclass of wavelet dictionaries. Sampling matrix 173 may include asampling operator determined by the transmitted waveform 152 associatedwith measurements 174 and the geometric arrangement of the independentsensors 184. Reconstructed image 160 may be transmitted by an electricalcommunication module 179.

FIG. 3, with reference to FIGS. 1 through 2, is a schematic diagramillustrating the options module 165, the estimation module 167, and thereconstruction module 169, according to an embodiment herein. In anembodiment, estimation module 167 includes first variable module 213configured to compute a first variable 211 based on a pre-selectednon-linear factor 209. In an embodiment, options module 165 determinesthe pre-selected non-linear factor 209. Options module 165 may furtherbe configured to determine a first pre-selected threshold 201 and asecond pre-selected threshold 203. In an embodiment, first variablemodule 213 is further configured to determine a multi-scale Gaussiantree structure 233 using a quad-tree decomposition of the electronicsource image model of the object/scene being sensed, the sampling matrix173, the pre-selected dictionary 161, and the measurements 174.

In an embodiment, the estimation module 167 may also include theparameter module 215 configured to determine a structural parameter 223using a parent-child relationship for each node in the tree structure233. The estimation module 217 may further include a loop module 217configured to control operation of the first variable module 213 and theparameter module 215 across the tree structure 233 until the structuralparameter 223 is lower than the first pre-selected threshold 201.Estimation module 167 may further include a second variable module 219configured to determine a second variable 234 using the first variable211, sampling matrix 173, variable selection operator 243, and thesecond pre-selected threshold 203. Estimation module 167 may alsoinclude a coefficient module 221 configured to determine coefficientvector 162 based on a Hadamard product of the first variable 211 and thesecond variable 234.

In an embodiment, estimation module 167 may be configured to (a)initialize x⁰=|h⁻¹(Ψ^(T)y)| and n=0, where h is the non-linear factor209, Ψε

^(m×d) is the sampling matrix 173, yε

^(m) are the measurements 174, and n is a loop counter; (b) calculatinga descent direction d^(n); (c) determining a step size λ; (d) computingx^(n+1)=x^(n)+λd^(n); (e) incrementing the loop counter n; (f) computingx* by repeating steps (b) through (e) until a norm of a steepest descentvector is smaller than the first pre-selected threshold 201; (g)computing z*=h(x*); (h) calculating

=diag(S_(Λ)|z*|) and Λ_(R)=diag(z*), where

$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix};} \right.}$(i) calculating u*=L⁻¹R, where L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is the noise factor 166; and (j) calculating c*=z*⊙u* where c*ε

^(n) is the coefficient vector 162.

In an embodiment, the reconstruction module 169 is further configured toreconstruct the image by calculating. I=Φc* where Iε

^(n) is the reconstructed image, and Φε

^(d×n) is the pre-selected dictionary 161.

FIG. 4, with reference to FIGS. 1 through 3, is a flow diagramillustrating a method 400 for reconstructing an image of an object,according to an embodiment. At step 402, method 400 furnishes, viasensing modality, measurements 174 of object 180 (modeled the global CGmodel) viewed by sensors 156. In an embodiment, sensors 156 areindependent from each other. In an embodiment, the image modeling module163 models the object 180 be imaged. At step 404, method 400 determinesthe sampling matrix 173, the pre-selected dictionary 161, and the noisefactor 166. In an embodiment, the sample matrix is determined using ageometric arrangement of sensors 156 collecting data of the object/scenebeing sensed. In an embodiment, option module 165 determines thesampling matrix 173, the pre-selected dictionary 161, and the noisefactor 166. In step 406, method 400 estimates the coefficient vector 162using the measurements 174, the sampling matrix 173, the pre-selecteddictionary 161, and the noise factor 166. In an embodiment, theestimation module 167 determines the coefficient vector 162. At step408, method 400 prepares a reconstructed image from the estimatedcoefficient vector 162 and the pre-selected dictionary 161 related tothe reconstructed image. In an embodiment, method 400 prepares thereconstructed image from the measurements 174, using the estimatedcoefficient vector 162 and the pre-selected dictionary 161 related tothe reconstructed image. In an embodiment, the reconstruction module 169prepares the reconstructed image.

FIG. 5, with reference to FIGS. 1 through 4, is a flowchart illustratinga method for computing the coefficient vector 162, according to anembodiment. At step 502, method 500 computes the first variable 211using the pre-selected non-linear factor 209, the multi-scale Gaussiantree structure 233 using a quad-tree decomposition of the electronicsource image model of the object/scene being sensed, the sampling matrix173, the pre-selected dictionary 161, and the measurements 174. In anembodiment the measurements 174 are in vector form. In an embodiment,the first variable module 213 computes the first variable 211 anddetermines the multi-scale Gaussian tree structure 233. At step 504,method 500 estimates the structural parameter 223 based on aparent-child relationship for each node in the tree structure. In anembodiment, the parameter module 215 determines the structural parameter223.

At step 506, method 500 repeats steps 502 and 504 across the treestructure until the parameter is lower than the first pre-selectedthreshold 201. In an embodiment, the loop module 217 controls operationof the first variable module 213 and the parameter module 215 across themulti-scale Gaussian tree structure 233 until the structural parameter223 is lower than the first pre-selected threshold 201. At step 508,method 500 computes the second variable 234 using the first variable211, the sampling matrix 173, a variable selection operator, and thesecond pre-selected threshold 203. In an embodiment, the second variablemodule 219 determines the second variable 234, the sampling matrix 173,the variable selection operator, and the second pre-selected threshold203. At step 510, method 500 computes the coefficient vector 162 usingthe Hadamard product of the first variable 211 and the second variable234. In an embodiment, the coefficient module 221 computes thecoefficient vector 162.

FIG. 6, with reference to FIGS. 1 through 5, is an embodiment forreconstructing the image of an object, according to an alternativeembodiment. At step 602 method 600 initializes x⁰=|h⁻¹(Ψ^(T)y)| and n=0where h is a nonlinear factor and Ψε

^(m×d) is the sampling matrix 173, yε

^(m) are the measurements 174, and n is a loop counter. At step 604,method 600 calculates a descent direction d^(n) and determines a stepsize λ. At step 606, method 600 computes x^(n+1)=x^(n)+λd^(n) andincrements the loop counter n. At step 608, method 600 computes x* byrepeating steps 604 and 606 until a norm of a steepest descent vector issmaller than a pre-selected threshold.

At step 608, method 600 may also compute z*=h(x*). At step 610, method600 calculates

=diag(S_(λ)[z*]) and Λ_(R)−=diag(z*) where

$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix}.} \right.}$At step 612, method 600 calculates u*=L⁻¹R where L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is the noise factor. At step 614, method 600 calculates c*=z*⊙u*where c*ε

^(n) is the coefficient vector. In an embodiment, estimation module 167performs steps 602 through 614 of method 600.

At step 616, method 600 calculates I=Φc* where Iε

^(n) is the reconstructed image and Φε

^(d×n) is the pre-selected dictionary 161. The pre-selected dictionary161 can include at least one class of wavelet dictionaries. The samplingmatrix 173 can include a sampling operator determined by a transmittedwaveform associated with the measurements 174 and a geometricarrangement of the independent sensors. In an embodiment, reconstructionmodule 169 reconstructs the image at step 616 of method 600.

FIG. 7, with reference to FIGS. 1 through 6, is a schematic diagramillustrating the multi-scale Gaussian tree structure 233 according to anembodiment. The multi-scale Gaussian tree structure 233 is based on aquad-tree decomposition of the image space. Each level of the quad-treedecomposition corresponds to a different scale of the waveletrepresentation of the image wherein each pixel at level i corresponds tod different wavelet coefficients corresponding to that location. Forexample, when choosing dictionary Φ to corresponding to orthogonalwavelets each node of the multi-scale Gaussian tree contains a vector ofsize d=4 corresponding to the familiar LL, LH, HL, HH sub-bands of thewavelet decomposition. In an embodiment herein, for simplicity, amulti-scale tree structure can be defined by the following equations:x(s)=A(s)x(par(s))+B(s)w(s)  (7)P _(x)(s)=A(s)P _(x)(par(s))A(s)^(T) +Q(s)  (8a)P _(x)(s,t)=Γ(s,sΛt)P _(x)(sΛt)Γ^(T)(t,sΛt)  (8b)where, s is node 11 of tree 13 and par(s) refers to the parent 15 ofnode 11;w(s)˜

(0,1),Q(s)=B(s)B^(T)(s); the symbol sΛt refers to the closest commonancestor to nodes s and t on a tree; and where Γ is a state transitionmatrix defined recursively as: Γ(s,s)=1,Γ(s,t)=A(s)Γ(par(s),t). Thesecoarse-to-fine dynamics implicitly define a probabilistic graphicalmodel in which edges between parent-child nodes represent a jointlyGaussian dependency structure.

The statistics of wavelet coefficients can be modeled effectively by theCG model and related observations in sea-clutter and radar images, asshown in Eq. (6). Given Eq. (6) and property (ii), the covariance of acoefficient and its parent is as follows:cov[c(s),c(par(s))]=

[h(x(s))h ^(T)(par(x(s)))]⊙cov[u(s)u ^(T)(par(s))]

The fact that the wavelet coefficients of images across scales tend tobe decorrelated can be enforced by constraining the Gaussian field u tobe corresponding white noise process across the multiscale tree.Furthermore, the variance of the wavelet coefficients across scalesfollows a self-similar structure. These properties can be accommodatedby modeling:u(s)=D(s)ç(s) ç(s)˜

(0,1) such that: D(s)=2^(−γm(s))where, m(s) is the scale corresponding to node s, and γ>0 is a constant.From Properties (i)-(iii) it follows that the variance of c iscontrolled by the u-field.

In spite of the above decorrelation across scales, the strongnon-Gaussian local dependence of wavelet coefficients that is observedin natural images is captured by the tree dependency structure in the x(and correspondingly, via (i), the z) random field. Finally as mentionedabove, the sparse structure of the wavelet coefficients is enforced bythe non-linearity h.

In using the Graphical CG model for the purposes of solving Eq. (1),simplifying assumptions can be made: (1) A(s)=A·I_(d) and B(s)=B·I_(d)∀swhere, A,Bε

, and I_(d) denotes the identity matrix of size d, and (2) the treestructure corresponding to x is homogeneous. From these assumptions, Eq.(7) can be written as:

$\begin{matrix}{{P_{x}(s)} = \frac{B^{2}}{1 - A^{2}}} & (9)\end{matrix}$Letting A≡μ and and B≡√{square root over (1−μ²)}, P_(x)(s)=I_(d)∀s.Thus, a single parameter μ controls the inter-scale dependency structureof the wavelet coefficients.

Though many different sparsity-inducing non-linearities h can beincorporated, in an embodiment herein, the following non-linearity canbe useful in terms of analytical properties such as smoothness withrespect to the sparsity controlling parameter α:h(x)=√{square root over (exp(x/α))}  (10)

Eq. (10) shows that the sparsity level of the signal is inverselyproportional to α. Thus two parameters, α and μ, control the propertiesof the Graphical CG distribution.

Given the CG model in Eq. (6), the pdf (probability density function) ofthe wavelet field c is given by:

$\begin{matrix}{{P(c)} = {\int_{\;}^{\;}{\frac{1}{\sqrt{2\pi}{\sum }z}{\exp\left( {{- \left( {c/z} \right)^{H}}{\sum^{- 1}\left( {c/z} \right)}} \right)}{{P_{z}(z)} \cdot {\mathbb{d}z}}}}} & (11)\end{matrix}$The structure of the CG distribution results from the summation of acontinuum of different Gaussians with different scales (variances)—eachof which is weighted by the profile P_(z)(z). Different choices of theweighting profile P_(z)(z) result in synthesis of different kinds pdfs(each with different kinds of heavy-tailed behavior)—including many ofthe well-known distributions in statistics, for example, but not limitedto, the Laplacian, Generalized Gaussian, and alpha-stable distributions.

A statistical model with prior distribution π(θ) is said to behierarchical Bayes if it can be decomposed in terms of conditionaldistributions

π(θ|θ₁), π(θ₁|θ₂), . . . , π(θ_(t−1)|θ_(t1)), and marginal π(θ_(t)):π(θ)=∫π(θ|θ₁)π(θ₁|θ₂) . . . π(θ_(t−1)|θ_(t1))π(θ_(t)).dθ ₁ . . . dθ_(t)  (12)The CG distribution inherently has a hierarchical structure in thatconditioning on the latent variable z Gaussianizes the wavelet field c.A fully Bayesian non-parametric approach to statistical inference undera hierarchical Bayesian prior model involves drawing samples from theposterior distribution:P(c,z|y)∝P(c,z)P(y|c)=P(c|z)P(z)P(y|c)  (13)via Markov Chain Monte Carlo (MCMC) and other sampling methods. Theintermediate estimation of the latent variable, the z-field, is referredto herein as Type-II estimation, whereas the estimation of the primaryparameter of interest c is referred to herein as Type-I estimation.

An embodiment focuses on a particular family of prior distributions;i.e., the Graphical CG model described above. The embodiment performsinference by a sequence of MAP estimates starting with the latentvariable:

-   -   1) Step 1: Perform Type-II MAP estimation problem:        z*=argmax _(z) log P(z|y)    -   2) Step 2: Perform Type-I MAP estimation problem conditioned on        z*:        c*=argmax _(c) log P(c|y,z*)

An embodiment solves the optimization problem presented in Eq. (4)wherein a Graphical CG model is placed on the unknown coefficients c.The first step is to estimate the non-Gaussian component z (Type-IIestimation), followed by the estimation of u (Type-I estimation). Withrespect to Type II estimation, given that z=h(x), it suffices toestimate the Multi-scale Gaussian Tree random process x. This can beperformed by recourse to the Bayesian-MAP strategy:x*=argmax _(x) log P(x|y)  (14)argmax _(x) log P(y|x)+log P(x)  (15)From Eq. (1) and Eq. (6):y 32 {tilde over (Ψ)}(h(x)⊙u)+v where, {tilde over (Ψ)}≡ΨΦ

y=A _(x) u+v  (16)where, A_(x)={tilde over (Ψ)}H_(x)H _(x) =diag(h(x))From Eq. (16):y|x˜

(y;0,A _(x) P _(u) A _(x) ^(T)+Σ_(v))x˜

(x;0,P _(x))where,

(w;μ,Σ) denotes a Gaussian with mean μ and covariance matrix Σ. Thus Eq.(16) is equivalent to:x*=argmax _(x) f(x)  (17)where:f(x)=y ^(T)(A _(x) P _(u) A _(x) ^(T)+Σ_(v))⁻¹ y+log det(A _(x) P _(u) A_(x) ^(T)+Σ_(v))+x ^(T) P _(x) ⁻¹ x  (18)

In an embodiment, two types of steepest descent approaches areconsidered: Gradient and Newton based steepest descent. The Gradientdescent involves a first-order Taylor series expansion of the costfunction Eq. (18) wherein only the Gradient vector of f has to becomputed in each iteration. The Newton descent approach, which involvesthe second order Taylor approximation of f, in addition entails thecalculation of the Hessian matrix. For Newton iterations, once theHessian is calculated via equation, the closest psd (positivesemi-definite) approximation to the Hessian is computed and used forcalculating the descent direction as described above. Given a matrix X,an approximation to its closest psd can be calculated as follows:M=QLQ ^(T)  (19)where,

-   L≡diag(max(real(Λ_(x)),0))-   Q≡√{square root over (Q_(X) ^(T)Q_(X))}-   Λ_(X)≡diagonal matrix of eigenvalues of X-   Q_(X)≡matrix of eigenvectors of X

Performing the closest psd approximation of the Hessian matrixX=∇²f(x^(n)), which effectively constrains the search direction to anelliptical ball defined by the matrix M, furnishes descent directionsthat can improve upon the gradient direction.

With respect to Type I estimation, let x* be the solution of Eq. (17)obtained via Type II estimation. Thus, the estimate of the non-Gaussiancomponent of the unknown coefficient c is z*=h(x*). The goal of theType-I estimation is to estimate the corresponding u field. TheBayesian-MAP strategy is used to generate the optimality criterion tosolve:

$\begin{matrix}{u^{*} = {\arg\;{\max_{u}{\log\;{P\left( {u❘y} \right)}}}}} & {{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}(20)} \\{= {{\arg\;{\max_{u}{\log\;{P\left( {y❘u} \right)}}}} + {\log\;{P(u)}}}} & {(21)} \\{= {{\underset{u}{\arg\;\min}\left( {y - {A_{x^{*}}u}} \right)^{T}{\sum\limits_{\upsilon}^{- 1}\;\left( {y - {A_{x^{*}}u}} \right)}} + {u^{T}P_{u}^{- 1}u}}} & {(22)}\end{matrix}$Eq. (22) reduces to solving the following equation with respect to u:Λ_(L)({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)})Λ_(R) u=Λ _(L){tildeover (Ψ)}^(T)Σ_(v) ⁻¹ y  (23)where,Λ_(L)=Λ_(R) =diag(z*)

The sparse structure of z* can allow pruning the number of equations tobe solved:

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) u=

{tilde over (Ψ)} ^(T)Σ_(v) ⁻¹ y  (24)

=diag(S_(λ) [z*])  (25)

$\begin{matrix}{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ \begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix} \right.} & (26)\end{matrix}$

S_(λ) is a variable selection operator with respect to threshold λ. Thethreshold λ is determined from the histogram of z*. The resulting Type-Iestimation methodology can be summarized as follows:

1. Let z* be the solution of the Type-II estimation procedure givenmeasurements y, sampling matrix Ψ, and dictionary Φ

2. Calculate

=diag(S_(λ)[z*]) for a threshold λ (in eq. (25)) and let Λ_(R)=diag(z*)

3. Solve for u:u*=L ⁻¹ Rwhere,L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R)R=

{tilde over (Ψ)} ^(T)Σ_(v) ⁻¹4. Return u*

The matrix

prunes out the rows corresponding to values of the z that are too small(as determined by threshold λ).

The HB-MAP methodology can be summarized as follows:

1. Initialize parameters

2. Perform Type-II estimation to obtain z*, the MAP estimate of the zfield

3. Given z* perform the approximate EM methodology to estimate theoptimum parameter μ

4. Iterate between steps (2)-(3) until convergence of μ

5. Perform Type-I estimation to obtain u*

6. Return the optimum estimate of the image: I*=

c*, where c*=z*⊙u*

In addition to employing the Type-I and Type-II estimation proceduresdescribed above, the HB-MAP methodology also employs an approximate EMmethodology that enables the refinement of parameter μ that controls thenon-Gaussian inter-scale dependencies between the wavelet coefficients.Let μ_(k) be the value of μ in the k^(th) iteration, and let x*_(k) bethe corresponding Type-II MAP estimate that solves (14). A 2^(nd) orderTaylor expansion of the cost function f about x*_(k) can be computed asfollows:f(x)≈f(x* _(k);μ_(k))+[∇f(x* _(k);μ_(k))]^(T)(x−x* _(k))+0.5(x−x* _(k))^(T){∇² f(x* _(k);μ_(k))}⁻¹(x−x* _(k))  (27)≈0.5(x−x* _(k))^(T){∇² f(x* _(k);μ_(k))}⁻¹(x−x* _(k))  (28)since ∇f(x*_(k);μ_(k))=0 by definition of the fact that x*_(k) is alocal optimum point, and f(x*_(k);μ_(k))≈0 by assumption. Thusp(x|y,μ_(k))≈

(x;x*_(k),∇²f(x*_(k);μ_(k))).

The Q-function of the EM methodology is:Q(μ;μ_(k))=−

_(p(x|y,μ) _(k) ₎[log p(x|y,μ)]  (29)Q(μ;μ_(k))=n/2 log(2π)+1/2 log det(P _(x)(μ))+1/2(x* _(k))^(T) P _(x) ⁻¹(μ)x* _(k)+1/2tr(P _(x) ⁻¹(μ)∇² f(x*_(k);μ_(k)))  (30)

and the gradient function of Q is as follows:∇_(μ) Q(μ;μ_(k))=−1/2(∇_(μ) P _(x)(μ))vec(P _(x) ^(−T)(μ))+1/2(∇_(μ) P _(x)(μ))vec((P_(x) ^(−T)(μ))x* _(k)(x* _(k))^(T) P _(x) ^(−T)(μ))+½(∇_(μ) P _(x)(μ))vec((P _(x) ^(−T)(μ))∇² f(x* _(k);μ_(k))P _(x)^(−T)(μ))  (31)The optimum parameter μ can be found as follows:μ*=argmin _(μ) Q(μ;μ_(k))  (32)where Eq. (32) can be solved by a simple gradient descent by using thegradient expression given in Eq. (31).

The computationally most intensive part of the methodology is thecalculation of the Hessian matrix and (to a lesser extent) the Gradientvector that is required in the steepest descent methodology (Type IIestimate). Direct implementation of the Hessian matrix and Gradientvector equations without exploiting the sparse structure of the matricescould be computationally expensive due to the presence of massive tensorproducts resulting from Kronceker operations on high dimensionalmatrices. The computational bottleneck can be overcome as follows. Thegradient vector is given by the following equation:∇f(x)=2P _(x) ⁻¹ x+G _(x)(Ψ^(T)

Ψ^(T))v _(x)  (33)where,G _(x)=∇_(x) H _(x)[{(P _(u) H _(x))

I _(n) }+{I _(n)

(P _(u) P _(x))}]  (34)v _(x) =vec((M _(x) ^(T))⁻¹)−[(M _(x) ⁻¹ y)

(M _(x) ⁻ y)]  (35)

$M_{x} = {{{A_{x}P_{u}A_{x}^{T}} + {\sum\limits_{\upsilon}\mspace{14mu}{{and}\mspace{14mu}{\nabla_{x}{H_{x}\left( {i,j} \right)}}}}} = \left\{ {{\begin{matrix}\frac{\partial{h\left( x_{i} \right)}}{\partial x_{i}} & {j = {{i\left( {n + 1} \right)} - n}} \\0 & {else}\end{matrix}P_{x}^{- T}} = {{P_{x}^{- 1}\mspace{14mu}{and}\mspace{14mu} M_{x}^{- T}} = M_{x}^{- 1}}} \right.}$

The following is a procedure for calculating the gradient in Eq. (33).

1. for i=1 to n

2. r=G_(x)(i,:)

3. index=find locations k where r(k)≠0

4. g=0

5. for j=1 to |index|

6. k=index(j), k_(mod)=mod(k,n), k_(s)=(k−k_(mod))/n

7. v₁={tilde over (Ψ)}(k_(s),:) v₂={tilde over (Ψ)}(k_(mod),:)

8. g=g+v_(x) ^(T)vec(v₂v₁ ^(T))

9. endfor

10. endfor

11. g=g+2P_(x) ⁻¹x

12. Return g

The Hessian form can be expressed as follows:∇² f(x)=2P _(x) ⁻¹ x+Σ _(i=1) ⁴ M _(i) [N ₁ +N ₂ ]+Q ₁ +Q ₂  (36)where,M ₁=∇² H _(x)[(P _(u) H _(x))

I _(n) ₂ ]({tilde over (Ψ)}^(T)

{tilde over (Ψ)}^(T)

I _(n))  (37)M ₂=∇² H _(x) [I _(n)

(P _(u) H _(x))

I _(n)]({tilde over (Ψ)}^(T)

{tilde over (Ψ)}^(T)

I _(n))  (38)M ₃ =∇H _(x) [I _(n)

(P _(u) ^(T)

{tilde over (K)} _(nn))][I _(n) ₂

(∇H _(x))^(T)]({tilde over (Ψ)}^(T)

{tilde over (Ψ)}^(T)

I _(n))  (39)M ₄ =∇H _(x)(I _(n)

P _(u) ^(T))({tilde over (K)} _(nn)

I _(n))[I _(n) ₂

(∇H _(x))^(T)]({tilde over (Ψ)}^(T)

{tilde over (Ψ)}^(T)

I _(n))  (40)furthermore,N ₁ =vec(M _(x) ⁻¹)

I _(n)  (41)N ₂=(M _(x) ⁻¹ y)

(M _(x) ⁻¹ y)

I _(n)  (42)

$\begin{matrix}{Q_{1} = {- {\begin{bmatrix}{\left\{ {{{\overset{\sim}{\Psi}}^{T}\left( {M_{x}^{- 1}y} \right)}\left( {M_{x}^{- 1}y} \right)^{T}\overset{\sim}{\Psi}} \right\} \otimes} \\\left( {{\overset{\sim}{\Psi}}^{T}M_{x}^{- 1}\overset{\sim}{\Psi}} \right)\end{bmatrix}}}} & (43)\end{matrix}$

$\begin{matrix}{Q_{2} = {\begin{bmatrix}{{K_{nn}\left\{ {\left( {{\overset{\sim}{\Psi}}^{T}M_{x}^{- 1}\overset{\sim}{\Psi}} \right) \otimes \left( {{{\overset{\sim}{\Psi}}^{T}\left( {M_{x}^{- 1}y} \right)}\left( {M_{x}^{- 1}y} \right)^{T}\overset{\sim}{\Psi}} \right)} \right\}} +} \\{\left( {{{\overset{\sim}{\Psi}}^{T}\left( {M_{x}^{- 1}y} \right)}\left( {M_{x}^{- 1}y} \right)^{T}\overset{\sim}{\Psi}} \right) \otimes \left( {{\overset{\sim}{\Psi}}^{T}M_{x}^{- 1}\overset{\sim}{\Psi}} \right)}\end{bmatrix}}} & (44)\end{matrix}$

=∇_(x) H _(x)[{(P _(u) H _(x))

I _(n) }+{I _(n)

(P _(u) H _(x))}]  (45)

A property that can enable efficient calculation of the above matrixproducts is to preserve sparsity of the intermediate matrices productscalculated from left to right. First, by construction the matrices∇H_(x) and ∇²H_(x) are sparse. Given this, the following Hessianprocedure can allow the sparsity of the intermediate matrix products tobe preserved while producing the desired matrix products (for example,M_(i)N_(j)∀i,j).

1. {A_(i)ε

^(n) ^(i) ^(xn) ^(i+1) }_(i=1) ^(k) are matrices to be multiplied; n₁=n

2. for i=1:n

3. r=A_(i)(i,:)

4. for j=1:k

5. r=prod(r,A_(j))

6. end

7. D(i,:)=devec(r)*vec(N)

8. end

9. Return D

where, N is either N₁ or N₂ (in Eq. (41-42)) and prod( ) refers tospecialized matrix products that use algorithmic techniques similar tothat in the gradient calculation algorithm.

In terms of flop count, the Gradient calculation takes O(n.S) and theHessian calculation takes O(n^(s).S) (where integer S depends on matrixP_(u); for a block-diagonal matrix P_(u), S is bounded by a smallconstant) on a single-core CPU (Central Processing Unit). Each of theiterations in the Gradient procedure and the Hessian procedure areindependent of each other and thus can be performed in any order. Thiskey property allows for parallel processing to be efficiently exploitedin multicore CPU or GPU (Graphical Processing Unit) platforms.

FIG. 8, with reference to FIGS. 1 through 7, is a graph illustratingperformance results for high-sparsity (α=0.2) according to anembodiment. Monte Carlo trials in which 16×16 images are reconstructedbased on measurements from a Radon transform at eight differentuniformly-spaced angles are conducted. Given the received signal, theembodiments herein are used to reconstruct the image which is thencompared with the original image. In all the simulations, Φ is adictionary of Biorthogonal 1.1 wavelets. Three different versions of theembodiments herein are shown: Newton HB-MAP 21, Gradient HB-MAP 23, andPerfect Knowledge HB-MAP 25. Perfect Knowledge HB-MAP refers to amethodology where perfect knowledge of the z-field is known, whichtherefore serves as a benchmark for the performance of othermethodologies.

Newton 21, Gradient Descent 23, and Perfect Knowledge 25 are performanceresults from execution of the embodiments herein. Line 23 corresponds tothe Gradient descent iteration and line 21 corresponds to the Newtoniteration. The simulation results show that the Gradient descentiteration and the Newton iterations both have the maximum MSE at 2trials. Line 25 illustrates the Perfect Knowledge HB-MAP 25corresponding to the case where the z-field is known exactly, as insimulated Monte-Carlo trials. Thus, the Perfect Knowledge HB-MAPprocedure only applies the Type-I estimation to estimate the image. Theperformance of Perfect Knowledge HB-MAP is therefore a lower bound whichno other methodology can surpass.

The threshold λ that is used in the Type-I estimation procedure isdetermined in these experiments by the coordinate of the histogram ofthe estimated z field below which 0.1 of the samples is contained. Line27 corresponding to Il-Is, line 29 corresponding to CoSaMP, line 31corresponding to ROMP, and line 33 corresponding to BCS illustrateperformance results from high sparsity-based convex optimizationmethodologies. Line 35 corresponding to Fourier illustrates performanceresults from a standard Fourier-based reconstruction. Simulated imageshave high-sparsity (i.e. α=0.2). The high sparsity case is more typicalof natural images. The performance of the embodiments herein exceedsthat of other options.

FIG. 9, with reference to FIGS. 1 through 8, is a graph illustratingsimulated performance for 16×16 images with low-sparsity (α=1)reconstructed based on measurements from a Radon transform at eightdifferent uniformly-spaced angles, according to an embodiment. Line 21Lcorresponding to the Newton low, line 23L corresponding to the GradientDescent low, and line 25L corresponding to the Perfect Knowledge lowperformance measurements are taken based on execution of the embodimentsherein. Il-Is low 27L, CoSaMP low 29L, ROMP low 31L, and BCS low 33L areperformance results from sparsity-based convex optimizationmethodologies. Line 35L corresponding to the Fourier low illustratesperformance results from a standard Fourier-based reconstruction.Simulation results illustrates that the Newton HB-MAP low and GradientHB-MAP low, according to the embodiments describes herein, are among thebest perform better than other methodologies. Newton HB-MAP low andGradient HB-MAP low give identical performance for low-sparsity case.

FIG. 10, with reference to FIGS. 1 through 9, illustrates an exampleimage processed according to embodiments herein. FIG. 10 illustrates theperformance of the Newton HB-MAP algorithm, Gradient HB-MAP, the FBP,and the various CS algorithms when applied to inverting the Radontransform, taken at 15 uniformly spaced angles of a 32×32 segment of theBarbara image—and in which measurement noise is added at eachobservation channel (sensor). The quality of the reconstruction has beenmeasured using the VSNR quality measure. Image 45 is reconstructed usingNewton HB-MAP. Image 51 is reconstructed using Fourier. Image 46 isreconstructed using Gradient HB-MAP. Image 47 is reconstructed usingIl-Is. Image 49 is reconstructed using BCS. Image 53 is reconstructedusing ROMP. Image 55 is reconstructed using CoSaMP. The images arederived when applied to inverting the Radon transform, taken at fifteenuniformly spaced angles, of a 32×32 segment 43 of Barbara image 41 andin which 60 dB of measurement noise is added at each observation channel(sensor). The quality of the reconstruction has been measured using theVSNR quality measure.

FIG. 10 illustrates that the Newton- and Gradient-HB-MAP algorithms aresignificantly better than all of the CS algorithms and especially thetraditional FBP approach to imaging.

FIG. 11, with reference to FIGS. 1 through 10, illustrates an exampleimage processed according to embodiments herein. FIG. 11 illustrates theperformance of the Newton- and Gradient-HB-MAP algorithms when appliedto inverting the Radon transform, taken at 15 uniformly spaced angles,of a 64×64 segment of the Barbara image—and, again, in which measurementnoise is added at each observation channel (sensor). Image 63 isreconstructed using Gradient HB-MAP methodology when applied toinverting the Radon transform, taken at fifteen uniformly spaced angles,of a 64×64 segment 61 of Barbara image 41 and in which 60 dB ofmeasurement noise is added at each observation channel (sensor). Image62 is reconstructed using Newton HB-MAP. Image 65 is reconstructed usingIl-Is. Image 67 is reconstructed using BCS. Image 69 is reconstructedusing Fourier. Image 71 is reconstructed using CoSaMP. Image 73 isreconstructed using ROMP. When compared to other algorithms, both theNewton- and Gradient-HB-MAP outperforms every other CS algorithm and isalso once again much better than the standard FBP algorithm. The resultsin FIG. 10 and FIG. 11 also demonstrate the robustness of the HB-MAPalgorithm in noise compared to the prior art.

In an exemplary embodiment, the various modules described herein andillustrated in the figures, for example systems and devices illustratedin FIGS. 1 through 3, are embodied as hardware-enabled modules and maybe configured as a plurality of overlapping or independent electroniccircuits, devices, and discrete elements packaged onto a circuit boardto provide data and signal processing functionality within aspecial-purpose computer. An example might be a comparator, inverter, orflip-flop, which could include a plurality of transistors and othersupporting devices and circuit elements. The modules that are configuredwith electronic circuits process computer logic instructions capable ofproviding digital and/or analog signals for performing various functionsas described herein. The various functions can further be embodied andphysically saved as any of data structures, data paths, data objects,data object models, object files, database components. For example, thedata objects could be configured as a digital packet of structured data.The data structures could be configured as any of an array, tuple, map,union, variant, set, graph, tree, node, and an object, which may bestored and retrieved by computer memory and may be managed byprocessors, compilers, and other computer hardware components. The datapaths can be configured as part of a computer CPU that performsoperations and calculations as instructed by the computer logicinstructions. The data paths could include digital electronic circuits,multipliers, registers, and buses capable of performing data processingoperations and arithmetic operations (e.g., Add, Subtract, etc.),bitwise logical operations (AND, OR, XOR, etc.), bit shift operations(e.g., arithmetic, logical, rotate, etc.), complex operations (e.g.,using single clock calculations, sequential calculations, iterativecalculations, etc.). The data objects may be configured as physicallocations in computer memory and can be a variable, a data structure, ora function. In the embodiments configured as relational databases (e.g.,such Oracle® relational databases), the data objects can be configuredas a table or column. Other configurations include specialized objects,distributed objects, object oriented programming objects, and semanticweb objects, for example. The data object models can be configured as anapplication programming interface for creating HyperText Markup Language(HTML) and Extensible Markup Language (XML) electronic documents. Themodels can be further configured as any of a tree, graph, container,list, map, queue, set, stack, and variations thereof. The data objectfiles are created by compilers and assemblers and contain generatedbinary code and data for a source file. The database components caninclude any of tables, indexes, views, stored procedures, and triggers.

Some components of the embodiments herein can include a computer programproduct configured to include a pre-configured set of instructionsstored in non-volatile memory, which when performed, can result inactions as stated in conjunction with the methods described above. Forexample components of FIG. 1 may include a computer program product. Inan embodiment, imaging device 150 may include a computer program. In anembodiment, pre-processor module 156 may include a computer program. Inan embodiment, options module 165 may include a computer program. In anembodiment, estimation module 167 may include a computer program. In anembodiment, reconstruction module 169 may include a computer program. Inan embodiment, first variable module 213 may include a computer program.In an embodiment, parameter module 215 may include a computer program.In an embodiment, loop module 217 may include a computer program. In anembodiment, second variable module 219 may include a computer program.In an embodiment, coefficient module 221 may include a computer program.In an example, the pre-configured set of instructions can be stored on atangible non-transitory computer readable medium or a program storagedevice. In an example, the tangible non-transitory computer readablemedium can be configured to include the set of instructions, which whenperformed by a device, can cause the device to perform acts similar tothe ones described here.

The embodiments herein may also include tangible and/or non-transitorycomputer-readable storage media for carrying or having computerexecutable instructions or data structures stored thereon. Suchnon-transitory computer readable storage media can be any availablemedia that can be accessed by a special purpose computer, including thefunctional design of any special purpose processor, module, or circuitas discussed above. By way of example, and not limitation, suchnon-transitory computer-readable media can include RAM, ROM, EEPROM,CD-ROM or other optical disk storage, magnetic disk storage or othermagnetic storage devices, or any other medium which can be used to carryor store desired program code means in the form of computer executableinstructions, data structures, or processor chip design. Wheninformation is transferred or provided over a network or anothercommunications connection (either hardwired, wireless, or combinationthereof) to a computer, the computer properly views the connection as acomputer-readable medium. Thus, any such connection is properly termed acomputer-readable medium. Combinations of the above should also beincluded within the scope of the computer-readable media.

Computer-executable instructions include, for example, instructions anddata which cause a special purpose computer or special purposeprocessing device to perform a certain function or group of functions.Computer-executable instructions also include program modules that areexecuted by computers in stand-alone or network environments. Generally,progran modules include routines, programs, components, data structures,objects, and the functions inherent in the design of special-purposeprocessors, etc. that perform particular tasks or implement particularabstract data types. Computer executable instructions, associated datastructures, and progran modules represent examples of the program codemeans for executing steps of the methods disclosed herein. Theparticular sequence of such executable instructions or associated datastructures represents examples of corresponding acts for implementingthe functions described in such steps.

The techniques provided by the embodiments herein may be implemented onan integrated circuit chip (not shown). The chip design is created in agraphical computer programming language, and stored in a computerstorage medium (such as a disk, tape physical hard drive, or virtualhard drive such as in a storage access network). If the designer doesnot fabricate chips or the photolithographic masks used to fabricatechips, the designer transmits the resulting design by physical means(e.g., by providing a copy of the storage medium storing the design) orelectronically (e.g., through the Internet) to such entities, directlyor indirectly. The stored design is then converted into the appropriateformat (e.g., GDSII) for the fabrication of photolithographic masks,which typically include multiple copies of the chip design in questionthat are to be formed on a wafer. The photolithographic masks areutilized to define areas of the wafer (and/or the layers thereon) to beetched or otherwise processed.

The resulting integrated circuit chips can be distributed by thefabricator in raw wafer form (that is, as a single wafer that hasmultiple unpackaged chips), as a bare die, or in a packaged form. In thelatter case, the chip is mounted in a single chip package (such as aplastic carrier, with leads that are affixed to a motherboard or otherhigher level carrier) or in a multichip package (such as a ceramiccarrier that has either or both surface interconnections or buriedinterconnections). In any case, the chip is then integrated with otherchips, discrete circuit elements, and/or other signal processing devicesas part of either (a) an intermediate product, such as a motherboard, or(b) an end product. The end product can be any product that includesintegrated circuit chips, ranging from toys and other low-endapplications to advanced computer products having a display, a keyboardor other input device, and a central processor, and may be configured,for example, as a kiosk.

The embodiments herein can include both hardware and software elements.The embodiments that are implemented in software include but are notlimited to, firmware, resident software, microcode, etc. Furthermore,the embodiments herein can take the form of a computer program productaccessible from a computer-usable or computer-readable medium providingprogram code for use by or in connection with a computer or anyinstruction execution system. For the purposes of this description, acomputer-usable or computer readable medium can be any apparatus thatcan comprise, store, communicate, propagate, or transport the programfor use by or in connection with the instruction execution system,apparatus, or device.

The medium can be an electronic, magnetic, optical, electromagnetic,infrared, or semiconductor system (or apparatus or device) or apropagation medium. Examples of a computer-readable medium include asemiconductor or solid state memory, magnetic tape, a removable computerdiskette, a random access memory (RAM), a read-only memory (ROM), arigid magnetic disk and an optical disk. Current examples of opticaldisks include compact disk-read only memory (CD-ROM), compactdisk-read/write (CD-R/W) and DVD.

A data processing system suitable for storing and/or executing programcode will include at least one processor coupled directly or indirectlyto memory elements through a system bus. The memory elements can includelocal memory employed during actual execution of the program code, bulkstorage, and cache memories which provide temporary storage of at leastsome program code in order to reduce the number of times code must beretrieved from bulk storage during execution.

Input/output (I/O) devices (including but not limited to keyboards,displays, pointing devices, etc.) can be coupled to the system eitherdirectly or through intervening I/O controllers. Network adapters mayalso be coupled to the system to enable the data processing system tobecome coupled to other data processing systems or remote printers orstorage devices through intervening private or public networks. Modems,cable modem and Ethernet cards are just a few of the currently availabletypes of network adapters.

A representative hardware environment for practicing the embodimentsherein is depicted in FIG. 12, with reference to FIGS. 1 through 11.This schematic drawing illustrates a hardware configuration of aninformation handling/computer system 1500 in accordance with anexemplary embodiment herein. The system 1500 comprises at least oneprocessor or central processing unit (CPU) 110. The CPUs 110 areinterconnected via system bus 112 to various devices such as a randomaccess memory (RAM) 114, read-only memory (ROM) 116, and an input/output(I/O) adapter 118. The I/O adapter 118 can connect to peripheraldevices, such as disk units 111 and storage drives 113, or other programstorage devices that are readable by the system. The system 1500 canread the inventive instructions on the program storage devices andfollow these instructions to execute the methodology of the embodimentsherein. The system 1500 further includes a user interface adapter 119that connects a keyboard 115, mouse 117, speaker 124, microphone 122,and/or other user interface devices such as a touch screen device (notshown) to the bus 112 to gather user input. Additionally, acommunication adapter 120 connects the bus 112 to a data processingnetwork 125, and a display adapter 121 connects the bus 112 to a displaydevice 123 which may be embodied as an output device such as a monitor,printer, or transmitter, for example. Further, a transceiver 126, asignal comparator 127, and a signal converter 128 may be connected withthe bus 112 for processing, transmission, receipt, comparison, andconversion of electric or electronic signals.

The foregoing description of the specific embodiments will so fullyreveal the general nature of the embodiments herein that others can, byapplying current knowledge, readily modify and/or adapt for variousapplications such specific embodiments without departing from thegeneric concept, and, therefore, such adaptations and modificationsshould and are intended to be comprehended within the meaning and rangeof equivalents of the disclosed embodiments. It is to be understood thatthe phraseology or terminology employed herein is for the purpose ofdescription and not of limitation. Therefore, while the embodimentsherein have been described in terms of preferred embodiments, thoseskilled in the art will recognize that the embodiments herein can bepracticed with modification within the spirit and scope of the appendedclaims.

What is claimed is:
 1. A method of reconstructing an image of an object,said method comprising: determining, by a plurality of sensors, awaveform based on said object, wherein said plurality of sensors viewsaid object from a plurality of directions; determining, by apre-processing module, a plurality of measurements of said object usingsaid waveform, wherein said plurality of measurements are arranged in avector form; determining, by an option module, a sampling matrix, adictionary, and a noise factor, wherein said sampling matrix representsa geometric arrangement of said plurality of sensors, and saiddictionary is pre-selected by said option module; estimating, by anestimation module, a coefficient vector using said measurements, saidsampling matrix, said noise factor and a global compound Gaussian prior;and reconstructing, by a reconstruction module, said image, using saidcoefficient vector and said dictionary.
 2. A method of reconstructing animage of an object, said method comprising: determining, by a pluralityof sensors, a waveform based on said object, wherein said plurality ofsensors view said object from a plurality of directions; determining, bya pre-processing module, a plurality of measurements of said objectusing said waveform, wherein said plurality of measurements are arrangedin a vector form; determining, by an option module, a sampling matrix, adictionary, and a noise factor, wherein said sampling matrix representsa geometric arrangement of said plurality of sensors, and saiddictionary is pre-selected by said option module; estimating, by anestimation module, a coefficient vector using said measurements, saidsampling matrix, and said noise factor; and reconstructing, by areconstruction module, said image, using said coefficient vector andsaid dictionary, wherein said estimating said coefficient vectorcomprises computing, by a first variable module: a first variable, usinga pre-selected non-linear factor; and a multi-scale Gaussian treestructure, using a quad-tree decomposition of said image, said samplingmatrix, said dictionary, and said measurements.
 3. The method of claim2, wherein said estimating said coefficient vectors further comprises:estimating, by a parameter module, a structural parameter based on aparent-child relationship for each node in said tree structure;repeating, by a loop module: said computing of said first variable andsaid multi-scale Gaussian tree structure, and said estimating of saidstructural parameter, across said tree structure, until said structuralparameter is lower than a first pre-selected threshold.
 4. The method ofclaim 3, wherein said estimating said coefficient vectors furthercomprises: computing, by a second variable module, a second variablebased on said first variable, said sampling matrix, a variable selectionoperator, and a second pre-selected threshold; and computing, by acoefficient module, said coefficient vector based on a Hadamard productof said first variable and said second variable.
 5. A method ofreconstructing an image of an object, said method comprising:determining, by a plurality of sensors, a waveform based on said object,wherein said plurality of sensors view said object from a plurality ofdirections; determining, by a pre-processing module, a plurality ofmeasurements of said object using said waveform, wherein said pluralityof measeaents are arranged in a vector form; determining, by an optionmodule, a sampling matrix, a dictionary, and a noise factor, whereinsaid sampling matrix represents a geometric arrangement of saidplurality of sensors, and said dictionary is pre-selected by said optionmodule; estimating, by an estimation module, a coefficient vector usingsaid measurements, said sampling matrix, and said noise factor: andreconstructing, by a reconstruction module, said image, using saidcoefficient vector and said dictionary, wherein said estimating saidcoefficient vector comprises: initializing, by said estimation module,x⁰=|h⁻¹(Ψ^(T)y)| and n=0, wherein x⁰ is an initial approximation of atemporary parameter, h is a nonlinear factor, Ψ is said sampling matrix,y is a vector comprising said measurements, and n is a loop counter;calculating, by said estimation module, a descent direction d^(n);determining, by said estimation module, a step size λ; and computing, bysaid estimation module, x^(n+1)=x^(n)+λd^(n), wherein x^(n) is an n^(th)approximation for said temporary variable and x^(n+1) is an (n+1)^(th)approximation for said temporary variable.
 6. The method of claim 5,further comprising: incrementing, by said estimation module, said loopcounter n; computing, by said estimation module, x*, by repeating saidcalculating said descent direction d^(n), said determining said stepsize λ, said computing x^(n+1), and said incrementing said loop countern until a norm of a steepest descent vector is smaller than apre-selected threshold, wherein x* is a temporary variable.
 7. Themethod of claim 6, further comprising: computing, by said estimationmodule, z*=h(x*); calculating, by said estimation module,

=diag(S_(λ)[z*]) and Λ_(R)=diag(z*), where$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix};} \right.}$ calculating, by said estimation module, u*=L⁻¹Rwhere L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is said noise factor; and calculating, by said estimation module,c*=z*⊙u* where c*ε

^(n) is said coefficient vector.
 8. The method of claim 7, wherein saidreconstructing said image further comprises: calculating, by saidreconstruction module, I=Φc* where Iε

^(n) is the reconstructed image and Φε

^(d×n) is said pre-selected dictionary.
 9. The method of claim 8,wherein said plurality of sensors are independent from each other,wherein said dictionary comprises at least one class of waveletdictionaries, and wherein said sampling matrix comprises a samplingoperator determined by a transmitted waveform associated with saidmeasurements and said geometric arrangement of said plurality ofsensors.
 10. An imaging device comprising: a plurality of sensorsconfigured to generate a waveform based on an object; a pre-processorconfigured to determine a plurality of measurements using said waveform,wherein said plurality of measurements are arranged in a vector form;and a central imaging device comprising: an option module configured todetermine a sampling matrix, a pre-selected dictionary, a pre-selectednon-linear factor, a first pre-selected threshold, a second pre-selectedthreshold, and a noise factor, wherein option module determines saidsampling matrix using a geometric arrangement of said plurality ofsensors; an estimation module configured to estimate a coefficientvector using said plurality of measurements, said sampling matrix, saidpre-selected dictionary and a global compound Gaussian prior; and areconstruction module configured to reconstruct image of said object,using said coefficient vector and said pre-selected dictionary.
 11. Animaging device comprising: a plurality of sensors configured to generatea waveform based on an object; a pre-processor configured to determine aplurality of measurements using said waveform, wherein said plurality ofmeasurements are arranged in a vector form; and a central imaging devicecomprising: an option module configured to determine a sampling matrix,a pre-selected dictionary, a pre-selected non-linear factor, a firstpre-selected threshold, a second pre-selected threshold, and a noisefactor, wherein option module determines said sampling matrix using ageometric arrangement of said plurality of sensors; an estimation moduleconfigured to estimate a coefficient vector using said plurality ofmeasurements, said sampling matrix, and said pre-selected dictionary;and a reconstruction module configured to reconstruct an image of saidobject, using said coefficient vector and said pre-selected dictionary,wherein said estimation module further comprises: a first variablemodule configured to: compute a first variable using a pre-selectednon-linear factor; and determine a multi-scale Gaussian tree structureusing a quad-tree decomposition of said sampling matrix, saidpre-selected dictionary, and said measurements.
 12. The device of claim11, wherein said estimation module further comprises: a parameter moduleconfigured to determine a structural parameter using a parent-childrelationship for each node in said multi-scale Gaussian tree structure;and a loop module configured to control operation of said first variablemodule and said parameter module across said multi-scale Gaussian treestructure until said structural parameter is lower than said firstpre-selected threshold.
 13. The device of claim 12, further comprising:a second variable module configured to determine a second variable usingsaid first variable, said sampling matrix, said variable selectionoperator and said second pre-selected threshold; and a coefficientmodule configured to determine said coefficient vector using a Hadamardproduct of said first variable and said second variable.
 14. An imagingdevice comprising: a plurality of sensors configured to generate awaveform based on an object; a pre-processor configured to determine aplurality of measurements using said waveform, wherein said plurality ofmeasurements are arranged in a vector form; and a central imaging devicecomprising: an option module configured to determine a sampling matrix,a pre-selected dictionary, a pre-selected non-linear factor, a firstpre-selected threshold, a second pre-selected threshold, and a noisefactor, wherein option module determines said sampling matrix using ageometric arrangement of said plurality of sensors; an estimation moduleconfigured to estimate a coefficient vector using said plurality ofmeasurements, said sampling matrix, and said pre-selected dictionary;and a reconstruction module configured to reconstruct an image of saidobject, using said coefficient vector and said pre-selected dictionary,wherein said estimation module is further configured to: initializex⁰=|h⁻¹(Ψ^(T)y)| and n=0, wherein x⁰ is an initial approximation of atemporary parameter, h is a nonlinear factor, Ψ is said sampling matrix,y is a vector comprising said measurements, and n is a loop counter;calculate a descent direction d^(n); determine a step size λ; andcompute x^(n+1)=x^(n)+λd^(n), wherein x^(n) is an n^(th) approximationfor said temporary variable and x^(n+1) is an (n+1)^(th) approximationfor said temporary variable.
 15. The device of claim 14, wherein saidestimation module is further configured to: increment said loop countern; and compute x* by repeating said calculating said descent directiond^(n), said determining said step size λ, said computing x^(n+1), andsaid incrementing said loop counter n until a norm of a steepest descentvector is smaller than a pre-selected threshold.
 16. The device of claim15, wherein said estimation module is further configured to: computez*=h(x*); calculate

=diag(S_(λ)[z*]) and Λ_(R)=diag(z*) where$\;{{S_{\lambda}\left\lbrack {z(i)} \right\rbrack} = \left\{ {\begin{matrix}1 & {{z(i)} > \lambda} \\0 & {{z(i)} \leq \lambda}\end{matrix};} \right.}$ calculate u*=L⁻¹R where L=

({tilde over (Ψ)}^(T)Σ_(v) ⁻¹{tilde over (Ψ)}+Λ_(L) ⁻¹Σ_(u) ⁻¹Λ_(L)⁻¹)Λ_(R) and R=

{tilde over (Ψ)}^(T)Σ_(v) ⁻¹ where vε

^(m) is said noise factor; and calculate c*=z*⊙u* where c*ε

^(n) is said coefficient vector.
 17. The device of claim 16, whereinsaid reconstruction module is further configured to reconstruct saidimage by calculating I=Φc* where Iε

^(n) is the reconstructed image and Φε

^(d×n) is said pre-selected dictionary.
 18. A non transitory programstorage device readable by computer, tangibly embodying a program ofinstructions executable by said computer to perform a method forreconstructing an image of an object, said method comprising:determining, by a plurality of sensors, a waveform based on said object,wherein said plurality of sensors view said object from a plurality ofdirections; determining, by a pre-processing module, a plurality ofmeasurements of said object using said waveform, wherein said pluralityof measurements are arranged in a vector form; determining, by an optionmodule, a sampling matrix, a dictionary, and a noise factor, whereinsaid sampling matrix represents a geometric arrangement of saidplurality of sensors, and said dictionary is pre-selected by said optionmodule; estimating, by an estimation module, a coefficient vector usingsaid measurements, said sampling matrix, and said noise factor; andreconstructing, by a reconstruction module, said image, using saidcoefficient vector and said dictionary, wherein said reconstructing saidimage further comprises using hierarchical Bayesian maximum a posterioriand using a global compound Gaussian model as a statistical prior.
 19. Anon-transitory program storage device readable by computer, tangiblyembodying a program of instructions executable by said computer toperform a method for reconstructing an image of an object, said methodcomprising: determining, by a plurality of sensors, a waveform based onsaid object, wherein said plurality of sensors view said object from aplurality of directions; determining, by a pre-processing module, aplurality of measurements of said object using said waveform, whereinsaid plurality of measurements are arranged in a vector form:determining, by an option module, a sampling matrix, a dictionary, and anoise factor, wherein said sampling matrix represents a geometricarrangement of said plurality of sensors, and said dictionary ispre-selected by said option module; estimating, by an estimation module,a coefficient vector using said measurements, said sampling matrix, andsaid noise factor; and reconstructing, by a reconstruction module, saidimage, using said coefficient vector and said dictionary, wherein saidestimating said coefficient vector comprises computing, by a firstvariable module: a first variable, using a pre-selected non-linearfactor; and a multi-scale Gaussian tree structure, using a quad-treedecomposition of said image, said sampling matrix, said dictionary, andsaid measurements.
 20. The program storage device of claim 19, whereinsaid estimating said coefficient vectors further comprises: estimating,by a parameter module, a structural parameter based on a parent-childrelationship for each node in said tree structure; repeating, by a loopmodule: said computing of said first variable and said multi-scaleGaussian tree structure, and said estimating of said structuralparameter, across said tree structure, until said structural parameteris lower than a first pre-selected threshold.
 21. The program storagedevice of claim 18, wherein said reconstructing said image furthercomprises using global compound Gaussian model as a statistical prior.